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Abstract 

This work analyzes the overall computational complexity of the stochastic Galerkin finite element method 
(SGFEM) for approximating the solution of parameterized elliptic partial differential equations with both 
affine and non-afHne random coefficients. To compute the fully discrete solution, such approaches employ 
a Galerkin projection in both the deterministic and stochastic domains, produced here by a combination 
of finite elements and a global orthogonal basis, defined on an isotopic total degree index set, respectively. 
To account for the sparsity of the resulting system, we present a rigorous cost analysis that considers 
the total number of coupled finite element systems that must be simultaneously solved in the SGFEM. 
However, to maintain sparsity as the coefficient becomes increasingly nonlinear in the parameterization, it 
is necessary to also approximate the coefficient by an additional orthogonal expansion. In this case we prove 
a rigorous complexity estimate for the number of floating point operations (FLOPs) required per matrix- 
vector multiplication of the coupled system. Based on such complexity estimates we also develop explicit cost 
bounds in terms of FLOPs to solve the stochastic Galerkin (SG) systems to a prescribed tolerance, which 
are used to compare with the minimal complexity estimates of a stochastic collocation finite element method 
(SCFEM), shown in our previous work [16]. Finally, computational evidence complements the theoretical 
estimates and supports our conclusion that, in the case that the coefficient is affine, the coupled SG system 
can be solved more efficiently than the decoupled SC systems. However, as the coefficient becomes more 
nonlinear, it becomes prohibitively expensive to obtain an approximation with the SGFEM. 

Keywords: stochastic Galerkin, stochastic collocation, sparse polynomial approximation, complexity 
analysis, explicit cost bounds, finite elements 


1. Introduction 

Nowadays, stochastic polynomial methods are widely used alternatives to Monte Carlo methods (see, 
e.g., [15]) for predicting the solution to physical and engineering problems described by parameterized partial 
differential equations (PDEs) with a finite number of random variables. In the last decade, two classes of 
such methods have been proposed that often feature much faster convergence rates: intrusive stochastic 
Galerkin (SG) methods and non-intrusive stochastic collocation (SC) methods. Both approaches typically 
employ a Galerkin projection in the physical domain, produced here by finite elements, and the resulting 
fully discrete approximations only differ in their choice of multivariate polynomials for the discretization in 
the stochastic domain. For details about the relations between these methods see [17, 18, 19, 21, 24], and 
for computational comparisons between the SG and SC methods see, e.g., [3, 13, 19]. 
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The focus of this paper is to provide explicit cost bounds for applying the stochastic Galerkin finite 
element method (SGFEM) to the solution of an elliptic PDE, with stochastic diffusion coefficient parame¬ 
terized by finitely many random variables. In particular, we focus on the cost of constructing isotropic total 
degree SG approximations when the coefficient has both affine and non-affine dependence on the parame¬ 
ters. Under very basic assumptions on the coefficient, the solution to this problem has been shown to have 
analytic regularity in the random variables (see [33]). As a result, SG approximations that employ a global 
orthogonal basis have been shown to be optimal projections in the L 2 sense, converging sub-exponentially 
with respect to the cardinality of the polynomial subspace [32]. However, the computational cost of solving 
the coupled SG system does not grow linearly in the cardinality of the given subspace. Therefore, the 
convergence estimates do not indicate the total complexity of obtaining the approximation for a prescribed 
tolerance. 

When the diffusion coefficient can be written as a sum of separable functions of the physical and random 
parameters, the coupled SG system can be written as a sum of Kronecker products of SG matrices and 
finite element stiffness matrices. For every SG matrix, each nonzero element leads to a nonzero block of the 
coupled SG system, where the size of the block equals the size of the finite element stiffness matrix. To solve 
the SG system, one must simultaneously solve all the coupled finite element problems. In the case that the 
coefficient is affine in the parameters, the number of nonzeros in each SG matrix is of order 0(M p ) [14], where 
M p is the cardinality of the isotropic total degree polynomial subspace of order pgN. Thus, a matrix-vector 
product involving the coupled SG system requires 0(JhM p ) floating point operations (FLOPs), where Jh 
is the number of physical degrees of freedom. Therefore, the work of solving the coupled SG system when 
employing an iterative method, e.g., conjugate gradient (CG), is of the order 0{JhM p Nf t ^ r ) where N^. 
is the number of iterations required to achieve a prescribed accuracy of the fully discrete approximation 
[3, 14, 34], 

On the other hand, when the diffusion coefficient is a general non-affine function of the random param¬ 
eters, the cost of obtaining an approximation with the SGFEM is not as obvious as before. In this setting 
we consider two cases, namely, the coefficient is: (1) a polynomial with respect to the random variables, 
and; (2) a transcendental function with respect to the random variables. In the first case, as we increase 
the order of the polynomial, the block-sparsity of the SG system decreases, resulting in a SG system that 
incrementally becomes block-dense [12, 14, 21, 23, 34, 35]. In the second case, a separable representation 
can be guaranteed with the use of an orthogonal expansion [36, 37], such that, substituting the expansion 
into the discretized PDE recovers the Kronecker product structure. However, when the expansion is not 
truncated, the SG system is known to be entirely block-dense [14, 23]. Without a priori knowledge on 
the exact sparsity of the SG matrices in this case, it was estimated that the complexity of matrix-vector 
multiplications of the SG Kronecker product system is between 0(Ji,M 2 ) and 0(JhM p ) [34]. As such, it is 
impossible to make a conclusive statement about the computational cost, and, more importantly, does not 
account for the two cases above, i.e., when the coefficient is possibly a truncated polynomial of fixed total 
degree r £ N such that 1 < r < oo. In these cases, the work of solving the coupled SG system with an 
iterative method is given by 0(JhM(p,r)N^ r ), where M(p,r) is the total number of 0{Jh) finite element 
problems that must be simultaneously solved. 

The key challenge of estimating the cost of solving the SG system when the coefficient is a (truncated) 
polynomial of finite order is to provide bounds on the block-sparsity of the matrix, i.e., nonzeros of the 
SG system. To achieve this, we provide a rigorous counting argument, which can be seen as a general¬ 
ization of results from [14], for the exact sparsity of the SG matrices for an arbitrary order orthogonal 
expansion of a non-affine coefficient. As a result, we are able to provide bounds for A A(p,r) of the order 
0(M p M r min{2 r , M|> }), where M r is the cardinality of the total degree polynomial subspace used in an 

orthogonal expansion of order r of the coefficient. This result provides sharper estimates than the bounds 
in the case of the full orthogonal expansion from [35] since it depends on the truncation order r, and allows 
us to estimate the total complexity of solving the coupled system for general non-affine coefficients. Since 
the counting argument for the sparsity of the SG system relies only on the SG discretization of an elliptic 
operator in terms of orthogonal polynomials, we note that this argument can be reused to estimate the 
complexity of solving similarly defined PDEs with this method. 

In addition, we also develop explicit cost bounds in terms of FLOPs to solve the SG system. Our 
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approach relies on e-complexity analysis, wherein we balance the errors arising from the approximation with 
the SGFEM and the iterative solver, e.g., CG, so as to ensure the solution to the fully discrete approximation 
achieves a given tolerance of e > 0. With this result, we are able to provide a direct comparison with e- 
complexity estimates for the stochastic collocation finite element method (SCFEM) in our previous work 
[16]. Finally, we present numerical results in agreement with the theoretical work estimates for both the 
SGFEM and SCFEM all cases described above. 

An outline of the paper is as follows. In §2, we provide a discussion on the model problem, and require¬ 
ments on the diffusion coefficient. In §3, we define the parameterized finite element and SG approximations, 
derive the SG system, and provide examples of the resulting linear systems that arise from the SG discretiza¬ 
tion with various coefficients. We then define the cost of solving the SG system and discuss preconditioning 
strategies. In §4, we derive the exact number of coupled finite element problems in the SG system and 
bounds on the sparsity in the non-afhne case, and present explicit cost bounds of the SGFEM. We also 
discuss the conditioning of the system in the non-affine case in order to provide a comparison with similar 
results from [29]. In §5, we briefly describe the SCFEM, and provide theoretical comparison with results 
from [16] in terms of minimum work to reach a given tolerance, both in the affine and non-affine cases. 
Finally, in §6, we present illustrative numerical examples corroborating our theoretical results. 

2. Problem setting 

We consider the simultaneous solution of the parameterized linear elliptic PDE: 

/ -V • (a(x,y)Vu(x,y)) = f(x) Vx G D, y £ T 

( u{x,y) = 0 Vx € dD, y € F ' ' 

where / £ L 2 (D) is a fixed function of x, D C R d , d = 1,2,3, is a bounded Lipschitz domain, and 
y(w) = (yi(w),..., vn{oj)) : fi —► r = Q R N is a random vector with u £ Q. and f l the set 

of outcomes. In this setting we assume the components of y have a joint probability density function 
g : r —>■ M + , with g{y) = J| i=1 Qi{Vi ) known directly through, e.g., truncations of correlated random fields 
[22] in (r, £>(r), g(y)dy), where B(y) denotes the Borel u-algebra on T and g(y)dy is the probability measure 
of y. We further assume that gi is an even weight function for each i = 1,... ,N. We require the following 
assumptions related to the continuity, coercivity, and holomorphic dependence of the coefficient a(x,y). 
Namely: 

(Al) There exist constants 0 < a m i„ < a max < 00 such that for all x £ D and y £T, a m i n < a(x, y) < a max . 

(A2) The complex continuation ofa(x , y), denoted a* : C N — > L°°, is a L°°(D)-valued holomorphic function 
on C N . 

The holomorphic dependence on y of the coefficient a(x, y) holds in many examples, including polyno¬ 
mial, exponential, and trigonometric functions of the variables yi,. ■■ ,yN shown below. 

Example 2.1 (The affine case). We consider an affine function of the random parameters, e.g., 

N 

a(x, y ) = a 0 (» + ykbk(x), x G D, y G T, (2) 

k =1 

where ao, {AfcjfeLi C L 2 (D) are such that a(x,y) satisfies (Al). Such examples include general Karhunen- 
Loeve expansions [22] or piecewise constant random fields. 

Example 2.2 (The non-affine, polynomial case). We consider a non-affine, polynomial function of the ran¬ 
dom parameters, e.g., 

a{x, y) = a 0 (x) + ^ y“c a (x), x£D,y£T, (3) 

l<|a:|<r 
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where a = (an,..., ajv) is a multi-index, |a| = a\ + • • • + a.^, y a = yf 1 ■ ■ ■ yffi, f < oo is the polynomial 
order ofa(x, y), and a o, {c a }\ a \<r C L 2 (D) are such that a(x, y) satisfies (Al). Examples include fixed-order 
Taylor or oHhogonal expansions of a general random field. 

Example 2.3 (The non-affine, transcendental case). We consider a non-affine, transcendental function of 
the random parameters, e.g., 


a(x,y) = a 0 (x) + g{x,y), x£D,y£T, (4) 

where ao , g C L 2 (D), and g(x , y) is a general transcendental function of x and y, such that a(x, y) satisfies 
(Al). Examples of g{x,y) include the sine, logarithm, or exponential functions of (2) or (3). 

Let L 2 ( T) be the space of square integrable functions with respect to the measure g{y)d.y and L“(r) be 
the space of essentially bounded functions, with the norm 

IMU~(r) := ess sup \u(y)\, 
yer 

where the essential suppremum is taken with respect to the weight g. By H~ X {D) we denote the dual 
of Hq(D), the space of square integrable functions in D having zero trace on the boundary and square 
integrable distributional derivatives. We will often use the abbreviation TL 2 to denote the space 

L 2 (T- Hq (D)) := : D x T —M : u strongly measurable and J \\u\\ 2 H i^g(y)dy < oo 

and 'H'ff to denote the space 


Lff (T; Hq (D)) :=\u:Dxr 


: u strongly measurable and ess sup ||u(-, y)\\H 1 (D) < 00 

yer 


For the space Hq(D) we have the energy norm |M|#i( D ) = ||Vv|| z , 2 (z>) , hence TL 2 is a Hilbert space with 
norm ||u||^ 2 = f r IMI^i^pdy. The stochastic weak form of problem (1) is given by: find u £ TL 2 such that 
\/v £ TL 2 


B[u,v]{y)g{y) dy = / F{v)g{y) dy, 


( 5 ) 


where 


B[u,v\(y) = [ , 
Jd 


a(x, y)Vu(x , y) ■ Vu( x, y)dx, 


F(v) = [ /( 

Jd 


x)v(x, y)dx. 


( 6 ) 


For convenience, we will often use the abbreviation B(y) = £>[•, -}(y) and suppress the dependence on x £ D 
in writing a{y) = a(-,y) and u(y) = u(-,y). It follows from (Al) that B{y) is a symmetric, uniformly 
coercive, and continuous bilinear operator on Hq(D), parameterized by y £ T, and B{y) induces the norm 


IMIB(y) := a(x, y)\X7u\ 2 dx. (7) 

Assumption (Al) and the Lax-Milgram lemma also ensure the existence and uniqueness of the solution u 
to (5) in H 2 . 

The convergence of the global stochastic polynomial methods used to approximate (1) exploits the 
uniform ellipticity of the coefficient a(y) and depends on the regularity of u(y) with respect to y. By Re(, 2 ) 
and hn(z) we denote the real and imaginary parts of z € C, and for 0 < <5 < a m ; n we define 

U(a,6) = {z£ C N : Re(o(a:,z)) > S,W x£ D}. (8) 
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If U (a, 8) ^ 0 for some 0 < S < a m i n , we say that a(x, z] is uniformly elliptic on the set U (a, 5) and we refer 
to Ufa, 5) as its domain of uniform ellipticity. For 7 = (71, , 7 ^) with 7 * > 1 Vi we denote the polyellipse 

£~t = (£) [zi G C : Re(zi) < 7 * + f l cos <j>, Im (zf) < 7i 7 * sin <f>, </> e [0, 27 r)|. 

In [33] it was shown that if a(y) satisfies (Al) and (A2), then for any 0 < S < a m j n there exists a 7 = 
(7i,...,7jv) with 7 , > 1 Vi such that C Ufa, 8). We can also similarly define the polydisc = 
£ C : \zi\ < 7 ,}, though, for arbitrary 0 < 8 < a m i n , it is not always possible to find a 7 with 
7 i > 1 Vi such that V-y C U(a,8). Figure 1 provides an illustration of this fact for various one-dimensional 
coefficients a(y), y £ C. Note that in the case of the 6 th degree polynomial and exponential random 
variables, no disc of radius 7 > 1 containing T = [—1,1] can fit in the region. The following theorem, proved 
in [33], shows the regularity of the solution u with respect to the parameterization. 




Figure 1: Domains of uniform ellipticity for some one-dimensional coefficients a(x,y) are indicated by the gray regions in each 
plot. The blue and red curves represent the maximal discs and ellipses, respectively, that can be contained in those domains, 
and the green lines represent the interval F = [—1,1]. 


Theorem 2.4. When the coefficient a(x,y) satisfies (Al) and (A2), so that for some 0 < 5 < a m [ n and 
7 = ( 71 ,..., 7 jv) with 74 > 1 V* we have £ C U(a, S), then the function z i —> u(z) from (1) is holomorphic 
in an open neighborhood of . 

This result states that a direct consequence of the uniform ellipticity of the function a{x , y) on the polyel¬ 
lipse £-y C C/(a, S) is that the solution u of (1) has analytically smooth dependence on the parameterization 
y. Theorem 2.4 is the key in motivating the construction of global stochastic Galerkin (SG) approximations 
to the solution u of ( 1 ), to be described in the following sections. 

3. Stochastic Galerkin finite element method 

In this section we define the SGFEM for constructing fully discrete approximations to the solution u of 
problem (1). This discretization employs mixed Galerkin projections in the spatial and parameter domains. 
In particular we rely on the finite element method for the spatial discretization, described in §3.1, and the 
stochastic Galerkin method for the parameter discretization, described in §3.2. In §3.3 we describe the linear 
systems that result from the SG discretization when Examples 2.1, 2.2, and 2.3 are used in problem ( 1 ). We 
then conclude in §3.4 with a discussion of the cost of solving the SG systems. 

3.1. Parameterized finite element approximation 

We briefly define the finite element method for obtaining a discretization of u from ( 1 ) over the spatial 
domain D. Let Th, be a triangulation of D with maximum mesh size h > 0, and VfJD) C Hq(D) a finite 


5 







































element space of piecewise continuous polynomials on Th parameterized by h —> 0. Let {4 , ji x )}'jLi denote a 
finite basis of Vh{D) of dimension J],. We can write the semi-discrete problem as: find Uh(y) G Vh{D) such 
that \/v G Vh(D) 


B[u h {y) 1 v)(y) = F(v), (9) 

where £>[-,-](y) and F(-) are defined in (6). For almost every y G T, problem (9) admits a unique solution 
of the form Uh(x, y) = J2j= i u j(y) < t > j( x )- We discretize problem (9) by defining, for i,j = 1,..., Jh , 

[A ] iJ {y)=B[4> j My), Fj = F{(j>i). (10) 

The coefficients u h (y) = \u\(y), u%{y), ... ,uj h (y)] T of Uh{x,y) are determined by solving the linear system 

Mv) u h(y) = F, (11) 

at fixed realizations of y G T. Here A (y) is symmetric and positive-definite so that (11) can be solved by 
iterative methods such as the conjugate gradient (CG) method. 


3.2. Stochastic Galerkin approximation with an orthogonal basis 

Based on the smoothness of the solution u to (1), characterized by Theorem 2.4, we now consider the 
construction of approximations to u in terms of global polynomials. Let A p C Nq be a finite set of multi¬ 
indices, e.g., having dimension #A p < oo, and define the space of polynomials Va p (T) = span {y u : u G A p }. 
A general global polynomial approximation problem can be framed in terms of solving for the ffA p stochastic 
degrees of freedom (SDOF) {tt p } p£ A p - When an interpolatory approach is used, the resulting systems of 
equations are decoupled finite element systems. When Galerkin projection with an orthogonal basis is used, 
the finite element systems are fully coupled, and must be solved simultaneously. Some isotropic examples 
of such index sets include 


A p = jpeN, 


Pi < , A™ = € Nq : Yj°n < P j , 

A™ = <> p G < : £ f{p n ) < f(p)) , f(p) = 


max 

l<i<N~ 


N 


n—1 


o, p = o 
i, p = i 
riog 2 (p)i, p> 2 


( 12 ) 


corresponding to the Tensor Products (TP), Total Degree (TD), and Smolyak (SM) polynomial spaces 
■P A tp(P), V\to(T), and ^^(P), respectively. When the solution u exhibits an anisotropic dependence on 
the parameters y, anisotropic weighted versions of the index sets defined in (12) can be introduced to further 
reduce the number of SDOF needed to approximate u at a desired accuracy [3, 26]. 


Remark 3.1 (Best M -term and quasi-optimal approximations). The optimal choice of A p would be the set 
A of cardinality Ad such that the corresponding approximation provides maximum accuracy out of all sets of 
size M. Such approximations are referred to as best AT-term approximations, and recent work has focussed 
on the construction of best M-term Taylor and Galerkin approximations [4, 6, T, 8, 33]. These approaches 
construct A by utilizing the largest M coefficients u p or sharp upper bounds of u p . However, in this effort 
we focus on analyzing the computational complexity of finding solutions to (5) in V\ (r) for a prescribed 
index set A p . 


For each n = 1,..., N, let {if Pn (2/n)}^ = i be a sequence of univariate polynomials over r„, orthonormal 

with respect to the L^ n (T n ) inner product. Then {'F p (y)}o<| p | with T p (y) := n^Li i’pniVn) is a sequence 
of multivariate polynomials over T, orthonormal with respect to the L^(r) inner product. In the case 
that g = \ for each n = 1,... ,1V, {Vw}p^=i anc i {^p}o<| P | are the univariate and multivariate Legendre 
polynomials, respectively. Given a specific choice of index set A p , it follows that {'P p } pe A p forms a basis of 


6 



V Ap (r) with dimension M p = dim('PA p (r)) = #A p . Hence, with {4>j}j^ ll as in §3.1 and {'kpjpgAp as above, 
we can now write the fully discrete stochastic Galerkin (SG) approximation as 


Jh 

tG, P (z,y)= 22 X] u i>p < M a 0'M y )’ 

p6A p j=1 


(13) 


whose coefficients can be found by solving the following coupled problem: find Uh, p £ 14(D) 0 Va (r) such 
that for all v £ Vh(D) 0 'PAp(r) 

E [B[u h>p ,v](y)\ = E [F(v)] , (14) 

where B[-, -](y) and F(-) are defined in ( 6 ). To form the linear system of equations resulting from the SG 
approximation given by (13), we let u h, P = [ui, p , ..., uj h , P ] T be the vector of nodal values of the finite 
element solution corresponding to the p-tli stochastic mode of Uh, P , and u/ ( p = [u/ liP ]p gAfi . Observe that 
when / is deterministic (H/pF;) = F,c>o,p for all i = 1 ,...,where <5o,p = 1 if p = 0 and <5o, p = 0 otherwise. 
Performing a Galerkin projection onto spanjd'pjpgAp for the solution of (14) yields the following system: 
for each p £ A p 


(^p(y)> A (y)^q(y)) u /i,<j = OMy), F ), (is) 

qeAp 

which can be written algebraically as a system of fully coupled finite element problems: for each p £ A p 

22 [K]p, q UA, q = F<5 0 p (16) 

qeAp 

with [K]p, q = (Tp(y), A(y)4' q (y)) and A (y) as given in (10). 

Remark 3.2. Typically matrix free methods are applied to solve (16) without ever explicitly forming K 
in memory, as described in [27]. When the resulting system is sparse, as a result of an affine coefficient 
a(x,y), e.g., Example 2.1, this can lead to computationally efficient solution strategies. However, these 
implementations rely on the fact that the coefficient a(x , y) can be written as a sum of separable functions 
of x and y, e.g., a(x, V) = E^Li b j(x)cj(y). For the transcendental function a{x,y) from Example 2.3, this 
may not be the case. Moreover, when K is block-dense, matrix-vector multiplications require approximately 
0(JhM p ) floating point operations (FLOPs), so that when iterative methods are used, the solution of the 
fully coupled finite element problems given in (16) becomes unfeasible. 

3.3. Representations of a(y) and the corresponding matrix K 

For a general coefficient a{x , y), the matrix K in (16) requires the storage of at most M p block matrices 
of the size and sparsity of A(y), i.e., O(JhM^) elements. However, in several specific cases the actual 
block-sparsity of K is much less. We recall the coefficient from Example 2.1, where K can be rewritten 

N 

[K]p,q = (4fp(y), 'kq(y))A q + 2_ y (yk^p{y ), 'bq(y))Afc, 

*=i 

with [Ao]ij = f D ao(x)Vtj>j(x) ■ V<fi(x)dx and [A*,]^- = f D bk(x)V<j>j(x) ■ V<fii(x)dx. If we let [Go] P ,q = 
('kp(y), 'Pq(y)) and [G fc ] Pi q = (y k ^p(y), ^q(y)), then K has a matrix representation, given by, 

N 

K = Go 0 Aq + 'y ( Gfc 0 Afc, (1^) 

fc =i 
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where A <g> B denotes the Kronecker product of A and B. We note that a similar construction can be 
optained for any coefficient a(x , y) which can be written as a sum of separable functions of x and y , such 
as the polynomial function of Example 2.2. 

However, when a(x, y) is not separable in x and y , this construction is no longer valid, and the resulting 
matrix K may be block-dense if we simply carry out the Galerkin projections and compute K directly. 
For certain special cases, e.g., when the diffusion coefficient is given by a log-transformed random field, 
the problem can be reformulated as a convection-diffusion problem, and the resulting system can be solved 
much more efficiently than the original problem [35]. In general, this reformulation is not applicable. 
Hence, for a general transcendental coefficient a(x,y), as given in Example 2.3, we project the coefficient 
onto an additional subspace r G No, in order to obtain a separable representation. To see this, 

define { v kr(y)}o<|r| to be the (infinite) basis of orthonormal polynomials of L 2 Q {T) as in §3.2. Then a(x,y) 
can be written as an expansion such that a(x,y) = X^o<|r| a r{x)^ r (y), where the coefficients a r (x) = 
(a(x,y), i S r (y)). Let A r be an index set of the type described in §3.2. Since infinite series representations 
are not practical in computations, we seek a truncation 

a r (x,y) := Y a r (x)^ r (y) (18) 

re A r 


in the subspace Vk t (r) for some r € No- When a r (x, y) ^ a(x, y), e.g., in the case that the projection order 
r is chosen to minimize error independent of the SG discretization, we let u r h denote the corresponding 
solution to the fully discrete SG approximation problem with a(x, y) replaced with a r {x, y). By substituting 
a r (x,y) into ( 6 ) we obtain 


/ f X Or(*)Ml/) ) ^4>j{ x ) ' ^M x )dx = Y [Arki’My), 

J D \ a / r6 A r 

[A r ]i,j = / a r (x)'V<j)j(x) ■ V<f>i(x)dx. 

Jd 


\reA r 


(19) 

( 20 ) 


Equation (19) represents an expansion of the stochastic finite element stiffness matrix A (y) and equation 
(20) represents the r-th mode of the expansion. Let u r hp = [u r i p ...., u r Jh p ] T denote the vector of nodal 
values of the finite element solution corresponding to the p-th stochastic mode of u r hp , and n r h p = [uj ( p ] pe A p . 
We substitute the expansion of A (y) into the Galerkin equations (15), to obtain the coupled system: for 
each p e A p 


X! X! [G’r]p,qA T .U ft) q — (TpjF), [G r ] p>q — (\l/p'I , q'I' r ). (21) 

A r qGAp 

Alternatively, similar to (17), we may define K r = ® A r to again obtain the coupled system of 

finite element problems: for all p £ A p 

Y \Kr]p, q <, q = F^o, P Vp G Ap. (22) 

qeA p 

Note that in forming K r , we now need only store the matrices {G r },.gA r and {A r } r£ A r , so that the efficient, 
matrix-free, solution strategies discussed in Remark 3.2 can be applied. 

Remark 3.3 (Projection and well-posedness). In the case that the coefficient is a transcendental function 
of the random variables, as in Example 2.3, there does not exist a r G No such that the projection (18) is 
exact. Due to the orthogonality of the basis, setting r = 2p in the construction of K r yields an entirely 
block-dense system [23] that is equivalent to (16), and computationally infeasible to solve. A more practical 
approach is to choose the expansion order 0 < r < 2 p, based on a-priori estimates of the error in the solution 
introduced by the truncation, so that the error when using the truncated expansion does not exceed that of 



the SG approximation. In this approach however, it becomes important to consider whether the truncated 
projection violates the well-posedness of (1) by failing to satisfy assumption (Al). One way to guarantee 
this is to choose f < r < 2p such that 


f := min{r <£ N 0 : ||a - a"||L~(r ; z,~(£>)) < a m in, W £ N 0 , v > r}. (23) 

An example of this problem can be seen in Figure 2 where for the function a(x, y ) = 0.1 + exp(2.5y), uniform 
ellipticity of the truncated projection a r (x,y) does not hold onF = [—1,1] for r < 4. 



Figure 2: Domains of uniform ellipticity for the total degree orthogonal expansions of order r of the one-dimensional coefficient 
a(x,y) = 0.1 + exp(2.5y), for y £ T = [—1,1] C R 1 , are indicated by the gray regions in each plot. The last plot shows the 
domain of uniform ellipticity of the original function a(x, y). The blue and red curves represent the maximal discs and ellipses, 
respectively, that can be contained in those domains, and the green lines represent the interval T. 


3.4 . Cost of solving the generalized SG system 

Without loss of generality, to solve the stochastic Galerkin system (22) for uj) , we use the precondtioned 
conjugate gradient (PCG) method, wherein, for the unpreconditioned CG method, we have the estimate 


*h,p 


— U 


r,(fc) | 
h,p 


K r ^ 


< 2 


- 1 


*~h,p 


— U 


r,(0)| 

h,p 


IK.,- 


(24) 


Here n r is the condition number of K r , u)jis the vector of the initial guess, and is the output of 

the fc-th iteration of the CG solver. The CG method is highly dependent on the conditioning of the system, 

r (k) 

and when n r is large, the number of iterations needed to reduce the error in u ;i ’ v p ; will also be significant. 
Hence we introduce the mean-based block-diagonal preconditioner (see, e.g., [27, 29]), 


P Go ® Aq, 


(25) 


with Aq and Go the matrices defined in (20) and (21) for r = 0 , respectively. 
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For r G A r , at every iteration of the CG method, or any iterative approach, each nonzero entry in 
each matrix G r implies a matrix-vector product of the form (\H p \t r g 1 I r T .)A T .pg fe \ where (’Fp’Fg’I'r) is a scalar 
quantity. Let nnz(A) denote the number of nonzeros of a matrix A, and define 

M(p,r) = ^ nnz(G r ) (26) 

re A r 


to be the total number of nonzeros in all of the matrices {G r } r£ A r at order p. With this in mind, an upper 
bound for the work in floating point operations (FLOPs) of solving (22) is given by 

W SG « 0(J h ) *M(p,r) * N'i L G , (27) 

where the term O(Jh) corresponds to the cost of a single finite element matrix-vector product, and _/V^ G is 
the number of iterations of the CG solver without a preconditioner. If we apply a preconditioner, in hopes 
to minimize A® G , we must also account for the added cost of applying the preconditioner at each iteration. 
With the mean-based preconditioner from (25), at each iteration we multiply an additional matrix of size 
JhM p x JhM p , but the matrix consists only of M p diagonal blocks. Here we assume that in finding the 
inverse of the preconditioning matrix P, a sparse approximation to P 1 is used, which we will denote P _1 . 
Such a decomposition can be found from, e.g., incomplete LU or incomplete Cholesky factorizations. Hence, 
for each iteration we require M p additional matrix-vector products of the size, and complexity, of the original 
finite element system, so the work estimate in FLOPs for the case of this preconditioner is given by 

W pSG « 0(J h ) * (M p + M ( P , r)) * TVPSG, ( 28 ) 

where A^® G is the number of iterations needed by the PCG method. Other preconditioners, such as the 
Kronecker product preconditioner suggested in [34] would require a different form of (28). 

Figure 3 displays the effect of fixing the projection order of the solution but increasing the order of 
the projection of the coefficient. In order to minimize the error of the projection, such a situation would 
be required if the coefficient is highly nonlinear and reflects the importance of considering Ai (p, r) in the 
computational cost of the SGFEM. 



Figure 3: Visualization of the number of nonzeros of a 165 x 165 SG matrix with elements [ K. r ]p,q = Sr*eA P r ]p >9 * J ^- r ‘ 
Each pixel represents a block finite element system when using a total degree projection of the solution of fixed degreee p = 3, 
and increasing the total degree of the projection of the coefficient, i.e., r = 0,1, 2, 3,4, 5. At r = 6, the matrix is entirely 
block-dense. 


4. Explicit cost bounds for the SGFEM 

The primary goal of this section is to estimate the algorithmic complexity required by the SGFEM to 
construct an approximation to (1) within a prescribed tolerance £ > 0. We assume a(x, y) is a general 
non-affine coefficient, as in Examples 2.2 and 2.3, satisfying assumptions (Al) and (A2). Let a r (x,y) be the 
orthogonal expansion of a(x, y), given by (18), of total degree r, i.e., a r (x, y) € 7 > A r (r) with A r = A™ from 
(12). We further assume that f < r < 2p, with r given in (23), so that a r (x,y) also satisfies (Al) and (A2). 
We will focus on the complexity of solving (22), when the stochastic discretization to (1) is performed in 
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V Ap (r) with A p = A™ from (12), i.e., in the space of total degree polynomials of order p, and the physical 
discretization is performed with the finite element method. These results are presented in the context of 
solving the linear system (22) with a PCG method when a zero initial vector is used to seed the solver. The 
results, however, can be generalized to other methods, such as preconditioned GMRES and other Krylov 
subspace methods. 

The results are organized as follows. In §4.1 we discuss the overall complexity of the matrix-vector 
products associated with solving (22) when using the SG matrix K r = A (g> A r . Our analysis 
extends the results of [14] in order to provide a bound on the block-sparsity of the SG system K r in the 
more general setting of a non-affine coefficient a r (x,y ), as given in Examples 2.2 and 2.3. In particular, 
for G r given in (21), we show that nnz(G r ) = 0(min{2l r 'l, M^ r y 2 ]}M p _^ r y 2 ]) for every r G A r , where 
M r = { N ^ r ) for r G No, when solving (22), so that the total complexity of the matrix-vector products with 
the Galerkin system is 0(JhM p M r min{2 r , Mr r / 2 i}). In §4.2, we perform an £-complexity analysis to derive 
the explicit cost bounds of the SGFEM using PCG, in terms of FLOPs as the tolerance £ —> 0. Finally, in 
§4.3 we discuss issues related to the conditioning of the SG system. 


4.1. Complexity of matrix-vector multiplications in the SG approximation 

In this section we provide rigorous bounds on the sparsity of the SG matrix K r from (22), for arbitrary 
0 < r < 2 p, and p G No- Our main result, given by Theorem 4.1, provides an exact count for nnz(G r ) in 
the general case |r| G No and N G N when the integrals [G r ] p , q = (4' T .4' p 4> q ) are defined in terms of even 
weight functions g. This result can be seen as an extension of estimates from [14], where bounds on the 
sparsity of G r were shown in the cases (i) |r| = 1 and N G N and (ii) N = 1 and r = r G No- We then 
provide upper bounds on the total number of nonzeros blocks A i(p,r) = YIvgA nnz (G r ) of the matrix K, 
from (22), both in the cases that a[x,y) is a finite order polynomial as in Example 2.2 and the case that 
a(x,y) is a transcendental function of the random variables, as in Example 2.3. Our first major result is 
summarized in the following Theorem: 


Theorem 4.1. Let A p and A r be the isotropic total degree index sets corresponding to the solution and 
coefficient, respectively, with p,r G No, and 0 < r < 2 p. If r G A r , and Qi are even for all i = 1,..., N, then 
for the matrix G r from (21) we have 


nnz(G r ) 



c(r,£) 


ffS(r,l) |r| even, I = |r|/2, 
2 #S(r,£) otherwise, 


(29) 


with S(r,£) = {sG N^ :\s\= s < r}, so that is equal to the coefficient oft e in the polynomial 

Mt) = ntiE;u P. Moreover, we have the following bound for nnz(G r ) ? i.e., 


nnz (G r ) < 2 min 



(GV+r|r|/2^ 


(N + P- \\r |/2l^ 


(30) 


so that 


A4 (p, r) < 2 ^ min 12 J , 
i=o ^ 


rr)} 


(N-l+j^ (N+p- 


(31) 


PROOF. For a given r G A r , we estimate the number of pairs (p, q) G A p x A p such that ('k r 4' p 4' q ) = 
nf=i<vv . tjj, Pt fqi ) yl 0. To do this, we extend the result of [14, Lemma 28] to a general matrix G r with 
|r* | G N. Since {4',.}, re A T . are orthonormal with respect to the even weight function p(y ) = IIi=i Pi(Ui)i we 
see that ('k T .4' p 4'q) ^ 0 if and only if (p, q) G O r , where 

Or = {(p, q) G A p x A p : | p t - q, \< n < pi + qi, 
and Pi + q, + ri is even Vi = 1,..., N}. 
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Therefore, to estimate the number of nonzeros in the matrix G r , we must estimate ^@ r . However, 0 r is 
different for each r £ A r . Even when r\ 1 r 2 £ A r are such that |r*i| = |r* 2 1, in general we do not have that 
#0 ri = #©r 2 - On the other hand, if r 2 is a permutation of rq, then it is easy to see that #O ri = #0 r2 
since A p is the isotropic total degree set. Also note that (’Pj.'Pp'I'q) = (d'r'I'qd'p) so if (p, q) £ 0 r then 
(q,p) £ O r as well. Note that 0 r can be rewritten 

Or = {(p, q) £ A p x A p : | p t - q,\ < n < pi + , 
and |Pi — qi\ + Ti is even Vz = 1,..., N}. 

Hence if (p, q ) £ 0 r , then we see that (p, q) must satisfy 

(i) | pi - q t | < n for all* = 1,..., N, 

(ii) fi <Pi+qi for all i = 1,..., N, and 

(iii) |pi — qi | + Ti is even for all z = 1,..., N. 

Note that when n = 0, we see that p,; = qi < p , and when > 0 we see that (i) and (iii) imply \pi — qi \ £ 
{0, 2,4,..., ri} for r.j even, and \pt — qi\ G {1, 3, 5,..., r.J for r,; odd. For each z = 1,..., N, let 
be the sequence defined by 

k (n) _ f 2?z + 1 ri odd, 

I 2 n ri even, 


so that fixing |p,; — <p| = k^ implies that conditions (i) and (iii) are met. 

To satisfy (ii) we must have r, ; < Pi + qi and to satisfy (i) and (iii) we must have \pi — q,;| = k^ n \ To 
avoid overcounting due to symmetry, we first fix possible values of pt and consider what qi must be. Let 
{4 7 °}k=o 2J be the sequence defined by 

, , r . . 7 » 

(n) _ r t + K i 


which we will refer to as the sequence of starting points for p,; corresponding to k[ n ^. Note that the starting 
points {si n) }£6 2j enumerate the integers between |"r*/2] and r*. Picking p 4 £ {s-"\s- n ^ + 1, ...,p} and 
Qi = Pi — k i we have 


Pi + qi = 2 pi - k^ > 2s l ( ” ) - /c[ n) 


A")' 


= 2 


- k\ n) = 


or Pi + qi > Ti, so that (ii) is satisfied. 

Since (i), (ii), and (iii) are satisfied by setting pi £ s^ + 1,... ,p} and qi = pt — k^ for a fixed 

0 < n < [z - i/2j, we count the number of admissible pairs for these choices. In N — 1 variables, the number 
of polynomials of total degree less than or equal to p — Pi is given by 

N - 1 + p - pi 
P~Pi 

where (?) = 0 if n < k or n, k < 0. To simplify notation, pick s, = s ^ (one of the starting points in the 
z-th direction) and fcj = k^ (its associated distance), where 0 < n < L r z/2J is fixed. To count the number 
of admissible pairs associated with the difference ki and starting point s^, we compute 


E 


Pi =s i 


nri-i'm 


N + p - sA 

p- Si ) 
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Define s £ with the .s, as above, then s corresponds to a possible combination of starting points in each 
direction. To estimate the number of polynomials associated with the starting point s, we compute 


p p—p i p—Pi - Pn-i 


EE- E 

Pi—Si P2—S2 Pw=s N 


'P-Pi - Pn\ 

:p - pi- pn ) 


fN+p- |s|\ 

1 p-\s\ )’ 


(32) 


where the sum easily follows by an induction argument and Pascal’s rule. 

Enumerating all of the pairs (p, q) £ Q r thus reduces to counting the number of possible combinations 
of starting points. Hence, in N dimensions we consider all such multi-indices of the {s ^ n> whose 
components sum to some integer |~|r|/2] < £ < |r|. For two multi-indicies s,r £ Nq , we say s < r if and 
only if Si < r, for all i = 1,..., N. Define the set S(r,£) = {s€ : |s| = £, s < r}, which corresponds 

to a particular slice of the desired set of starting points. To estimate #S(r,£), we consider the familiar 
counting argument of placing N bars among t stars with the added restriction that the numer of stars in the 
*-th bin not exceed r,. Such a problem can be reframed in terms of finding the coefficient c(r,£) of t e in the 
generating function P r (t) = JlE Y0j=■ Combining (32) and summing over £ between f|r*|/2”| < £ < \r\ 
we arrive at (29), where the coefficients c(r,£) = #S(r,£) when |r| is even and £ = |r|/2 (in this case the 
roles of pi and can not be reversed) and c(r ,£) = 2 #S(r,£) otherwise. 

Noting that £y^^ r ^ 2 -^S{ r ,£) is a change of coordinates of a total degree index set of order f|r*|/2”| 
intersected with the hyperrectangle {s £ : s < r} yields the bound 


M M 

E c ( r ^)< 2 E 

l=\\r\/2\ t=\\r\/2~\ 


#S(r,£)< 2 


^+r|r|/2^ 


On the other hand, from the generating function P r (t) we see that c(r,£) is bounded by (^) when |r| is even 
and £= |r|/2 and 2(y) otherwise. This follows from the fact that when k is the multi-index having |r| ones 
and the rest zeros, since £ < |r|, we have that #S(r,£) < #T(fc, £) where T(k, £) = {s £ Nq : |s| = £, s < k} 
and #T(k,£) is given by the coefficient of t e in P fe (i) = (l + t)l T 'l = from the binomial theorem. 

Then 

M kl 

E c ( r ,^) < 2 E c ( r ^) = 2 |r|+l , 

^lkl/21 e=o 

so that 


nnz(G T .) 


kl 

E 

fefkl/21 



N + p 
P~ 


£ 


< 2 min 



(N+ \\r\/2]^ 


^N + p- \\r\/2]^ 


showing (30). Substituting (30) into (26) shows the bound of M(p,r) from (31). □ 

We note that the bound of A i(p, r ) from (31) is an overestimate due to the particular form of (29), which 
is different for each r £ A r . As a consequence, we see that nnz(G r ) = 0(mm{2 ' r ', M^ r i mi }M p _^\ r \ m ) for 
r £ A r . For large N and small r, 2 r is smaller than M |>/ 2 ], however, as r —> oo the bound M|- r / 2 ] is sharper. 


13 




For the e-complexity analysis in the next section, we note that 


Mir,r) = E “" Z < G -) £ E Sminhn, (" + (J. r|/21 ) } ( N + P 

reA r re A r ' ' / J \ 

f N + \j/2] \ \ fN — 1 + j\ fN + p — fi/21 


re A r 


■-rM/2i 

N 


= 2 min < 2 J , . ... , , 

^ 1 ’ v N \ N-l A N 

i=o ^ 

MN j SN + p^^fN-l+j 


< 2 min < 2 


j=0 


iV- 1 


= 2 min < 2 


+ fr/2]\ 1 /7V+p\ /JV + f 
N ( { N { N 


(33) 


Figure 4 plots how sharply A i(p,r) is bounded by (31) and (33). We are also able to show that Theorem 
4.1 yields a sharp result in the case |r| = 1. 


Upper bounds on A4(p, r) with 


Upper bounds on <\4(p,r) with 



Figure 4: For r = p with p ranging from 0,1,..., 6 we plot for AT = 4 (left) and TV = 8 (right) the actual sparsity A4(p, r) 
given by (26) of the Galerkin system K r from (22) (blue), the bound on the sparsity from (31) (green), and the bound on the 
sparsity from (33) (red). 


Corollary 4.2. Under the same conditions in Theorem 4-T when r £ A r is such that |r| = 1, we have 



Corollary 4.2 is the result of [14, Lemma 28], and follows from the application of the exact formula for 
nnz(G T .) from (29). Here |r| = [|r|/2] = 1 is odd and S(r,l) has only one element S(r,l) = {s £ Nq : 
I s ! = 1, s < r} = {r}. Hence c(r, 1) = 2 #S(r, 1) = 2, and (34) is shown. We are also able to show that the 
formula for nnz(G r ) from (29) yields a result that is sharp in the case N = 1. 

Corollary 4.3. Under the same conditions in Theorem 4-1, when N = 1 and r = r £ No, we have 

(a) in case r = 2k, k £ No, 

{ (p — r + l)(r + 1) + k 2 , 0 < r < p, 

(p—k + 1) 2 , p+l<r<2p, (35) 

0, r > 2 p. 
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(b) in case r = 2k + 1, fc € No, 


! {p — r + 1 )(r + 1) + A: 2 + k, 0 < r < p, 

(p — k + 1 )(p — k), p + 1 < r < 2p, (36) 

0, r > 2p. 


Corollary 4.3 is the result of [14, Lemma 25], and its proof using Theorem 4.1 is included in the Appendix. 
In the remarks that follow, we make a distinction between the cases that a(x,y) is a polynomial of fixed 
degree F < oo, e.g., the coefficients from Examples 2.1 and 2.2, and that a(x , y) is a transcendental function 
of the random variables, e.g., the coefficient from Example 2.3. 


Remark 4.4. (Complexity of matrix-vector products for polynomial coefficients, see e.g., Examples 2.1 and 
2.2) From Corollary 4-2 and the work estimate (28) when using (25) as a preconditioner, we see that for 
coefficients that are affine functions of the random variables, e.g., Example 2.1, the complexity of a single 
PCG iteration is of the order 0(Jh{2M p + 27VM p _ 1 )) = 0(JhM p ), where M p = #Aj D = ( N ^ p ). On the 
other hand, when the coefficient a(x,y) is a polynomial function of the random variables, e.g., Example 2.2, 
having fixed order r € N, F < oo, we use Theorem f.l to obtain a different estimate. Since {4' r } T . e A F is a basis 
for the space there exits coefficients {a r (a;)} re A- such that a(x,y) = a r (x,y) = Ylre A - a r{x)^ r (y). 

With this representation, it is clear to see that substituting a(x,y) into (19) yields from (22), and 
Kr = K from (16). However, it is not clear how many of the coefficients a r (x) are identically zero. In 
this case, we can provide an upper bound on the block-sparsity of under the assumption that a r (x) ^ 0 
Vr £ A ¥ . Using the bound of (33), the complexity of a single matrix-vector product o/K f is of the order 
0(JhM p Mr min{2 r , M|>/ 2 -|}). Thus, when f is fixed, min{2 r , M^/z \} is a constant, and this estimate 
has the same asymptotic complexity as 0(JhM p ). 


Remark 4.5. (Complexity of matrix-vector products in the trascendental case, see e.g., Example 2.3) We 
recall the discussion of [34, Section 3-4]- There, the complexity of matrix-vector products with the SG 
system was estimated when a full orthogonal expansion is substituted into the SG discretization. This case 
corresponds to fixing the expansion order r = 2p following Remark 3.3. Assuming that nnz(G r ) = 0(M p ) or 
0(M p ), it was estimated in [35] that the cost of matrix-vector products involving K r is be between 0(JhM p ) 
and 0(JhM p ). However, the use of Theorem 4-1 allows us to consider the complexity in the case of truncating 

the expansion, where a sharper estimate can be obtained. Let T r := rifc=[>/ 2 l+i ^fr/2] = (^[v/^)’ 

which is bounded independent ofr, i.e., 


(N+\r/2] +1 

l fr/21+1 


fr/21+l 




as r —► oo, 


so that M r = T r M|- r / 2 ] < e N Mv r / 2 ft From (33), we see that A4(p, r) is of the order 0[M p M r M^ r /W) as 
p,r —> oo, since min{2 r , M|>/ 2 ]} —» M\ r /s\ as r 00 . When r = 2 p, this implies the complexity of matrix- 
vector multiplications involving K r is of the order 0(JhM p M r M^ r /2]) = 0(JhM p T r M ^ r/2 ^) = 0(JhM p ). 
On the other hand, when r = p, we see that the complexity of matrix-vector products with K r is order 
0{J h M p M r M ]rm ) = 0(J h T?Mf p/2} ) = 0(J h Mf pm ). 


4-2. e-complexity analysis of the SGFEM 

An estimate of the total complexity to obtain a fully discrete approximation of tolerance e > 0 with the 
SGFEM and PCG solver can be shown in four steps: 

1. Estimate the maximum mesh size h max and minimum polynomial order p m i n necessary in the finite 
element and SG discretizations, respectively, 

2. If projection of the coefficient is necessary, estimate the minimum projection order r m ; n , otherwise set 
rmin = r where F < 00 is the order of the coefficient, 

3. Estimate the minimum number of iterations k m \ n needed by the PCG solver, 
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4. Substitute h m ax , p m in> r m ; n , and fc m in into the cost (28) for h, p , r, and A^® r G , respectively. 

We proceed to estimate these parameters as follows. Denote by u r the corresponding solution of (1) when 
a r (x,y) is substituted in place of a(x,y ), and let u r hp be the approximation to u r hp found by PCG. Then 
the total error for the SG approximation satisfies the following bound: 

\\ u ~ KpIIh 2 < ll u — KK + IK — U h\\n 2 + || u h ~ KpIIh 2 + || u h,p ~ K,p 11-^2 • (37) 

SG(I) SG(II) SG(III) SG(IV) 

In this setting SG(I) is the approximation error using a truncated expansion of a(x,y), SG(II) is the 
discretization error induced by the finite element method, SG(III) is the SG error coming from the orthogonal 
expansion (13), and SG(IV) is the solver error resulting from the PCG method. We note that when the 
projection of the coefficient is exact, as discussed in Remark 3.3, the approximation error SG(I) is no longer 
present and u r hp = Uh, p - 

We start with bounding SG(III). Without loss of generality, it is reasonable to assume that since u r 
has a holomorphic dependence on z £ C N in an open neighborhood of the polyellips £-y from Theorem 2.4, 
then u r h does as well. Then, the following result, whose proof is found in [32], and immediately follows from 
classical spectral convergence results [9, 31], describes the convergence rate of the fully discrete solutions 
obtained by the SG method using a total degree approximation in 7 , a td(P): 

Proposition 4.6 (Convergence rate for the SG method). If Theorem 2-4 holds for the solution u r h to (9) 
with coefficient a r {x,y), and u r hp is the solution to (14) with A p the order p total degree index set, then 

IK - u h, P \\Hf> < Ci exp(—C 2 p) Vp € N, 

for some constants C\,C '2 > 0 independent of p. 

To investigate the error in SG(I), we note that since a(x,y) satisfies assumption (A2), the projection 
error in ’P A td(P) can be similarly estimated as 

11° - ; l~(d)) < C 3 exp(—C 4 r) Vr G N, (38) 

for some constants C 3 , C 4 > 0 independent of r. Hence, Vr G N, 

I \u-u r \\ n < ^^\\a-ar\\ Ll(T , L oo (D)) < M|^i C3 exp(-G 4 r) (39) 

^min ^min 

providing a bound for SG(I). For a bound of SG(II), we present the following convergence result reguarding 
solutions to the parameterized finite element problem, whose proof can be found in a number of standard 
texts on the theory of finite element methods, e.g., [2, 20]: 

Lemma 4.7. Let Th be a uniform finite element mesh over D with Jh = 0{h~ d ) degrees of freedom and 
h > 0. For the elliptic PDE (1) and y G T, when u r (y) G Hq(D) n H S+1 (D), the error from the finite 
element approximation is bounded by 

IK(y) - u r h {y)\\ H i {D ) < C FEM h s , 

where the constant C FEM > 0 is independent of h and y. 

For the treatment of SG(IV), we begin by defining B r (y) to be the corresponding bilinear operator 
in (6) with a(x,y) replaced with a r (x,y). Since both B(y) and B r (y) are symmetric, uniformly coercive 
and continuous bilinear operators on Hq(D), there exist a,/3 > 0 independent of y such that for every 
u,v £ Hq(D) 


\B r [u,v\(y) \ = 


a r (x, y)Vu ■ A/vdx 


id 


Ih 0 i ( d) < J^a r (x,y)\\7u\ 2 dx = 


< a ll u ll.ffi(D)IMIiJi(.D)) an d 


levy)’ 
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and similarly for B r (y) with the same a,/3, e.g., taking a to be the maximum and /3 to be the minimum 
in each case. Recall n r h = [u{ p ,... ,u r Jh p ] T , the vector of nodal values of the finite element solution 
corresponding to the p-th stochastic mode of u r hp , and u r h = [u Ap ] pgA . Then we have the following 
estimates expressing 


Continuity: IKJIk, = IKJIeib-^)] - II u hAw e ’ and ( 40 ) 

Ellipticity: sffi ||ufc iP || w? < |K, p |j E[Br(H)] = i( u h,p|| Kl . > ( 41 ) 


where ||u||^ = (u) r K r u is the K r matrix norm, and ||u||E[B r (p)] is the expectation of the energy norm (7). 
Given Proposition 4.6, Lemma 4.7, and the estimates from (39), (40), and (41), we can now provide the 
minimal projection orders p,r £ N for the SG approximation (13) and the coefficient (18), respectively, the 
maximum mesh size h for finite element method, and the minimum number of PCG iterations k necessary 
to ensure that the error in the SGFEM solution u r hp is less than the tolerance e > 0. 


Lemma 4.8. Let u £ L 2 g (T; Hq(D) n H S+1 (D)) be the solution to (1), u r hp be the solution to (14) with the 
coefficient a r (x, y), and u r h be the approximation of u r h found by PCG with a zero initial guess. Then, for 
e > 0, to ensure that ||u — u r h p ||^2 < £ we must choose h < h m ax , r > r m i n , p > p m in, and k > k m i n , where: 


bin ax — 


. dGpEM 
Pmin = log 


4Cf\^ 

£ 


= log 


4 C 5 \ ^ 


kmin — 


><* (®d) 


with C FEM > 0 the constant from Lemma f.7, Ci, C2, C3, C4 > 0 the constants from Proposition f.6 and 
(38), and, for a,/3> 0 from (40) and (41) 


C S = C 3 




mi 


with K r is the condition number of P 1 K r with P the mean-based preconditioner from (25). 


PROOF. Without loss of generality, we seek to bound the quantities SG(I)-SG(IV) from (37) each by e/4. For 
the error SG(I) we recall estimate (39) and solve for r. From Lemma 4.7, when u £ L 2 (F; Hg(D)nH s+1 (D)) 
we have that ||u r — u r h ||^2 < C FEM h s Vft > 0, and from Proposition 4.6 we have that \\u r h — u r hp \\-H™ < 
C i exp(— C 2 P) Vp G N, so that solving for h and p gives the desired maximum mesh size h m ax and minimum 
polynomial order p m ; n to bound SG(II) and SG(III) by e/4. Let u' hjl and u/'^ f be the coefficients of the 
exact SG solution u r h and the approximate SG solution u r h after k PCG iterations, respectively. Then 
from (24) and (41) we see that 


\^h,p nfop ||^2 — 


V/?' 


r,(fc)| 

l h,p u h,p I 


2 ( frO-r 1 \ 11 r r >(0)l 

TfVTrrr) l|u -‘-- u ^ IIk ' 


where is the initial guess used in CG and h r = cond(P x K r ) with mean based preconditioner P from 


(25). If we use the zero vector as the initial iteration in PCG, we have from (40) 


2 / frfrf — 1 


a 


ik r — 1 


IKp <,p\\ H 2 < [^- r + l ) IKpIk. < 2^ [frjff+i) bhAui- (42) 

Solving for k gives the minimum number of iterations fc m ; n required to ensure SG(IV) is bounded by e/4. □ 
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Given the necessary parameters from Lemma 4.8 to achieve ||w — ^hpWn 2 < £ i and the estimates on 
the computational complexity of one iteration in the PCG method from §4.1, we provide a bound on the 
minimal number of FLOPs required by the SGFEM when approximating (1). We split these results into 
the cases that the stochastic coefficient a(x,y) from (1) is: 

(i) an affine function of the random parameters, e.g., a(x,y) £ P,\. r . (P) with r = 1, as in Example 2.1, 

(ii) a non-affine polynomial of the random parameters, e.g., a(x,y) £ V (P) for some 1 < r < oo, as in 
Example 2.2, 

(iii) a non-affine, transcendental function of the random parameters, e.g., ci(x,y) ^ VA r (P) for any r £ N, 
as in Example 2.3, so that r must be chosen to satisfy r > r m i n from Lemma 4.8. 

The results are summarized in Theorems 4.9 and 4.10 next. 


Theorem 4.9. Let u £ L 2 g (F; Hq (D)nH s+1 (D)) be the solution to (1), and r be the smallest natural number 
such that a(x,y) £ Pa^(T). When f = 1, the minimum work (28) of solving (22) with PCG to a tolerance 
e > 0 can be bounded by 


W pSG < C 7 


3 c v 


26^(1 +TV) 


^1 + log 





(43) 


and when r > 1, the minimum work (28) of solving (22) with PCG to a tolerance e > 0 can be bounded by 


W pSG < C 7 


3G f 


Mw 1 ) 

+s(ifS) 


„N 


1 + fog 


3C'i \ °2? 
£ 


N 


i=o 


vrw< + r 




r , , i-i 


[(y)”’1 

)] 


(44) 


with C fem ,Ci,C 2 ,Cq as in Lemma /.8, C 7 > 0 independent of e, and h the condition number of the 
preconditioned system P _1 Kr = P _1 K, using the mean-based preconditioner from (25). 


PROOF. When a(x, y) £ 7^(1") we do not need to consider SG(I) from (37), and bound SG(II), SG(III), and 
SG(IV) by e/3. Hence, to minimize the error of the SG discretization, we choose p > p m ; n = log[(3Gi/e) 1 / c ' 2 ] 
which differs from the p m in stated in Lemma 4.8. For a uniform triangulation Tk, Jh = 0{h~ d ) so that 


Jh 




C 7 



(45) 


for some constant C 7 > 0 depending on the connectivity of the finite element mesh, but independent of e. 
In the case that r = 1, we substitute p rn m into (29) for the matrices G r having 0 < |r| < 1, and apply 
Stirling’s approximation to obtain 


-^Pmin “b -M(Pmin> 1) — 2 


N + Pmin 

TV 


2 TV 


N + Pmin ~ 1 

TV 


< 2e N (1 + TV) 1 + log 


3 Gi \ & 

£ 


N 


Similarly, when r > 1 we use the bound from (31) and Stirling’s approximation to obtain 

N 


fo^Pmin + M (Pmin, r) < e 


N 


1 + log 


3Ci \ 

£ 


+2 min ^ 2 J , 

i=o 


TV- 


+ ri/2i\ wjv-i+j 


tv / ( v tv — i 




r 1 i 


[(?)“] 

)] 
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As in the proof of Theorem 4.9, we substitute J/, max for J;, from (45), k m in for N^ r , and the bound for 
M pmlii + M.(p m i„, r m i n ) with p m ; n and A m i n from Lemma 4.8 into the cost (28) to complete the proof. □ 

Given Theorems 4.9 and 4.10 we see that the work of obtaining the fully discrete approximation using 
the SGFEM, with PCG as a solver, is asymptotically given by: 



(SG.l) (SG.2) (SG.3) 


where g(N) = N and k r = k if a{x,y) is an affine or non-affine, polynomial function of the random 
parameters of fixed order r < oo, e.g., Examples 2.1 and 2.2, and g(N) = 3 N when a(x,y) is a non-affine, 
transcendental function of the random parameters, e.g., Example 2.3, requiring a total degree orthogonal 
expansion of order r > r m i„ depending on e. Here, (SG.l), (SG.2), and (SG.3) correspond to the work 
required by the finite element, SG, and PCG methods, respectively. In particular, (SG.2) corresponds to 
the estimates for the sparsity of the Galerkin system K r from (22), and represents the number of coupled 
finite element systems that must be solved simultaneously by the PCG method. However, due to the bound 
(31), the asymptotic complexity in the cases that a(x,y) is affine or polynomial in y are the same. This 
does not imply that there is no need to consider the work estimates in these cases separately. Indeed, if 
a(x,y) is a polynomial having the representation ]CreA- a r (aO’I'r (u) where a r (x) ^ 0 for all r £ A ¥ , then 
the complexity of matrix-vector multiplications with K r is of the order 0( JhM v M T min{2 r , M \r/i\ })■ Here, 
the constant min{2 r , Af [r/al} grows rapidly with r, suggesting that higher order polynomial functions of 
y require additional cost. 
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4-3. Conditioning of the generalized SG system 

In this section, we discuss issues related to the conditioning of the linear system that results from the 
SGFEM discretization. We first recall [14, Theorem 10]: the eigenvalues of the matrices {G r } r£ A r from 
(21) lie in the interval [£ r ,S r ], where 

4- := min{4' r (y) : y € E r := max{^' T .(y) : y G (48) 

gm(i) j g a t ensor product grid of Gauss-Legendre quadrature points having m(l) = (m(?i),... ,m(l n )) points 
in each direction, and 1 is such that m(l n ) := p + [ fc,1 2 +1 ], n = 1,... ,1V. Since a r (x,y) satisfies (Al), the 
analysis of [29, Theorem 3.8] shows that the eigenvalues for the preconditioned system P -1 K r lie in the 
interval [1 — r r , 1 + T r \ where 

T.r = —— ^2 ^r||a r (a;)|| L oo (D) , f r = ^ H r ||a T .(x)|| i oo (£ , ) . (49) 

° min reAr a min reAr 

|- r |^0 M =^0 

As a result of (49), we see that in the case that the projection order r of the coefficient a r (x,y) depends 
on e, the condition number of the preconditioned system P -1 K r does as well through the number of terms 
in T r and f r . This should come as no surprise since even in the case of the Karhunen-Loeve expansion, the 
condition number of P _1 K depends on the number of terms in the truncated Karhunen-Loeve expansion 
which is chosen a-priori to minimize the error. 

5. Comparison with the SCFEM 

In this section we compare our explicit cost bounds for the SGFEM with the complexity estimates for 
the SCFEM developed in [16], when solving (1). The basic idea behind the SCFEM is to construct a fully 
discrete approximation in a subspace of 14(D) (g> L 2 g (T) by collocating senri-discrete solutions Uh from (9) 
on a deterministic set of points to obtain solutions {uh{-,yk)}^=i € 14(D). 

5.1. A generalized SCFEM using Lagrange interpolation 

To construct the stochastic collocation (SC) approximation, we consider a class of multi-index sets 
defined in terms of increasing functions m : —> N4 and g : —> N+. By m we specify the multivariate 

function m( 1) := (mi(Zi),--- ,mjv(ijv)) where each m n : N + N + is an increasing function, possibly 
different for each n = 1,..., N. Here the m n are referred to as growth functions, specifying how the number 
of points grows in the direction n. Associated with m n we define the left-inverse rrv n : N+ —► N+ by 
mn(q) = min{/c £ N+ : m. n {k) > g}, and let m) (q) = (rn\(qi),... ,m' N (qN))- In this case, we note that 
m' n = k and m n (m' n (k)) > k for each k € N + and n = 1,..., N. Given m and g we can define the 
multi-index set 

A r ,S = {<? e N + : g( m H<l + 1)) < L} , (50) 

to be used in constructing polynomial approximations. In particular, setting m n {j) = j for all j £ N+ and 
n = 1,... ,N, and defining 


N 

9tp (p) = max p„, STd(p) = V'bn - 1), 

l<n<N z ' 

n—1 


N 

9SM {P) = ^2 f(Pn), 

n—1 


(51) 


where f(p) is given in (12), and using the definition of A)”’ 9 from (50), we obtain the TP, TD, and SM 
index sets Aj p , A™, and A| M , respectively, given in (12). 
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We introduce a sequence of one-dimensional Lagrange interpolation operators : C'°(r„) - 

'Pm 7l {i n )-\iXn)- Then for v G C ,0 (T) the generalized multi-dimensional approximation operator X™’ 9 
C°(r) —> P A -, g (T) is given by 


/ N 

g(l)<L ie{0,l} w \n=l 

Construction of the approximation T™' 9 [v\(y) requires the independent evaluation of samples v(y) on a 
deterministic set of distinct collocation points G™’ 9 having cardinality Ml = ffG™' 9 - Applying X™' 9 \] from 
(52) to the semi-discrete solution Uh{x,y) of problem (9), we obtain the fully discrete SC approximation 

Uh,L(x,y) = I^' 9 [u h }(x,y). (53) 

One-dimensional abscissas. In this effort, we use three examples for constructing the fully discrete approx¬ 
imation. The first is that of a fully-nested rule constructed on the Clenshaw-Curtis choice of abscissas [5] 
with function 5td(p) from (51) and an isotropic growth rule m = (m,..., to) with to given by 

to(1) = 1, m{l n ) = 2 ln ~ 1 + 1 for l n > 1, (54) 

This is the classical Smolyak sparse-tensorization construction [30], and here the choice of m corresponds 
to a doubling growth rule that leads to a nested sequence of multi-dimensional grids, e.g., G™’ 9 C 
On the other hand, we can construct a sparse-Smolyak approximation on the Gauss-Legendre abscissas 
corresponding to the zeros of the Legendre polynomials {'F p }, as defined in §3. When the points are grown 
isotropically according to the linear growth rule with m = (to, ..., m) and to defined as 

m(l n ) = l n for l n G N, (55) 

and 5td(p) from (51), we obtain a grid that is not nested. Another construction that yields a sequence 
of nested grids is that based on the Leja points, defined as the sequence of points satisfying yu+i := 
argmax yerj , n U \ y — yj\ (see [10]). Here we take the Leja sequence of points with <?td from (51) and the 
isotropic linear growth function m from (55). 

5.2. Cost of solving the SCFEM systems 

To construct the fully discrete approximation with the SCFEM, we must solve M L distinct decoupled 
finite element systems, each dependent on a realization of the parameters yu G Q™' 9 for k = 1,..., Ml- 
Similar to the SGFEM, we can apply the PCG method to the solution of each system. Let iVj te (. be the 
number of iterations required by the CG method to solve the finite element system corresponding to yu 
and be the corresponding number of iterations when a preconditioner is used. We are interested in 

choosing a suitable preconditioning strategy to decrease the total number of iterations = Y^k =i 

required to obtain the fully discrete approximation Uh,L- We present a preconditioning strategy of choosing 

P 0 :=A(yi), (56) 

with A (y) from (10) the finite element stiffness matrix corresponding to the sample point y\ G G^' 9 , as the 
preconditioner for all of the individual finite element solutions. We refer to this choice of preconditioner as 
the level-zero preconditioner since it corresponds to the SC approximation at level L = 0. 

Since we apply CG to the solution of each individual finite element system, the work in floating point 
operations (FLOPs) required to obtain a fully discrete approximation with the SCFEM without a precon¬ 
ditioner is given by 


l?' 9 [v]{y)= E E 



w sc 


M l 

0 ( 4 )^^ 


k =1 


(57) 
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On the other hand, the “level-zero” preconditioner induces an an additional matrix-vector product requiring 
O(Jh) FLOPs per iteration when a sparse factorization of Po is used. Hence the work of solving (5) with 
PCG is given by 


WV sc *2*0(J h )*J2 N ^- ( 5g ) 

k =1 

Here, the reduction in work due to preconditioning will be seen in the number of iterations saved in each 

p (k) 

individual count N contributing to the sum. 


5.3. Comparing the explicit cost bounds of the SGFEM and SCFEM 

Given a particular “sparse” index set A p , we can find increasing functions m : Nf —> and g : —► 

N+, and L £ N such that A p = A™’ 9 from (50). In this setting, we can either use Galerkin projection or 
construct an interpolant to obtain an approximation to u in 7 7 A p (r < ). Let u Ap denote the Galerkin projection 
of u onto the space V Ap (T). Then we have the estimate 


ll M “ M A p ||z,2(r ; Hi(D)) 


vGH'(D)®V Ap (r) 


where C a > 0 depends on the coefficient a(x,y) and the bounds from assumption (Al). This estimate 
expresses optimality in the L 2 g {T) error of the Galerkin projection since C a does not grow with A p , and 
suggests that the Galerkin method is the best choice for approximating u in the space Va p (C). We can also 
define an interpolation operator I™’ 9 : C^P) —>• V A (T), and then we have the estimate 


\\u-T ?' 9 


< (Ca l min^ (r) ||«-n|| L oo (r;Hol(i ,)) 

= (C Al + l)||u - WA p ||z,oo(r ; jji(D)) 


(59) 


where C Al is the Lebesgue constant of Iff’ 9 . A good interpolant will be one for which C Al grows moderately 
with ffAff' 9 . For example, it is known (see [11, 16]) that for a one-dimensional Lagrange interpolation 
operator using a Clenshaw-Curtis rule, the Lebesgue constant is bounded by 2 log(m — 1) + 1, where m 
is the number of points. For the SC method, we define SDOF to be the total number of points needed 
to construct the approximation. From (59), if we only consider the number of SDOF needed to represent 
the solution, we expect the error for the Galerkin approximation to be much lower than the error in the 
interpolant. Indeed, this is reflected in our numerical results in Figures 7 and 9, and has been observed in 
previous comparisons [3, 13]. 

However, if we are willing to change the space A p , e.g., adding more interpolation points to gain a more 
stable interpolant by changing m or changing which points are included in the set Aff’ 9 by changing g , 
it might be possible to obtain an approximation with lower complexity to reach a given tolerance, despite 
having to solve more systems. Therefore, to properly compare the work involved in constructing u Ap and 
Iff’ 9 [u], we consider the computational complexity of both methods, not in terms of SDOF, but in terms of 
floating point operations (FLOPs). For a chosen A. pi this reduces to studying the complexity of the system 
resulting from Galerkin projections and the stability properties of the interpolant Iff’ 9 . 

Let Uh,L denote the numerical solution to the fully discrete approximation uh,l obtained with the SCFEM 
from (53) found by the PCG, and observe that we have a similar splitting to (37) for the error in the 
approximation 


||u - Uh,L \\ H 2 < ||U - ^11^2 + \\u h - U } i , L \\ U 2 + \\Uh,L - Uh,L \\ n 2 (60) 

sc(i) V scqi) / scqii) 

Note that unlike in the case of the SGFEM, the SCFEM does not require a further projection of the coefficient 
a(x,y), so that we do not need to consider the error ||it — u r ||^2 from (37). In addition, we do not need to 
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worry about well-posedness of the truncation as discussed in Remark 3.3. Similar to the complexity analysis 
for the SGFEM, we must choose h < /i max and L > L m j n so that the errors ||u — Uh\\n 2 Q from the finite 
element discretization and \\uh — Uh,L^n 2 e from the SC interpolation are both bounded by e/3. From this, 
a minimum tolerance r m ; n for the PCG solver can be derived and the maximum number of PCG iterations, 
with a zero initial guess, can be estimated [16]. In what follows, we present a result, whose proof can be 
found in [16, Theorem 4.7] that bounds the number of PCG iterations in the context of the work estimate 
(58). Using this estimate we can compare the cost in FLOPs for the SCFEM with the SGFEM results from 
Theorems 4.9 and 4.10 in the previous section. 

Theorem 5.1. Let u € L^(T; Hq(D) C\ H 3+1 (D)) be the solution to (1). Then for e > 0 arbitrary, the work 
of finding Uh,L, the approximation to the fully discrete SC solution Uh,L from (53) found by PCG, denoted 
by IU pSC , can be bounded by 

W*° < (^) 1 c 8 [log (2f2)]" [c. + ± log log (2fs)] ** (.1, 

* { log {¥) +c ' 1 + 2Jvloglog [k log ( 3 t £ )] } ■ 

Here Cfem from Lemma 4-7, CV from Theorem 4.9, and C$, Cg, Cio, Cu, Csc, o-nd r from [16, Theorem 4-7] 
are positive constants independent of e. Moreover, we define R = sup^gp n(y) where n{y) is the condition 
number of the preconditioned system Eg A (y) with Po from (56). 

Theorem 5.1 follows from the fact that -/V P ® C = ero , where N zero is the number of 

iterations needed by the SCFEM with a PCG and a zero vector initial guess. Substituting the bound on 
N zelo shown in [16, Theorem 4.7] into the work estimate (58) and using J/,. max = Cr(3C FEM /e) d ^ s as in 
Theorem 4.9, puts the result in terms of FLOPs. Given Theorem 5.1 we see that the work of obtaining the 
fully discrete approximation with the SCFEM with the PCG method is asymptotically bounded by: 



(SC.l) (SC.2) (SC.3) 


where (SC.l), (SC.2), and (SC.3) correspond to the work required by the finite element, SC interpolant, and 
PCG methods, respectively. Since the costs associated with the finite element discretization are the same 
for both methods, we focus only on the costs associated with the SC projection and the SC interpolation, 
coupled with the costs of the PCG method. In particular, as the work required by the SC approximation 
from (SG.2) of (47) has different bounds depending on whether the coefficient is a fixed order polynomial or 
is given by a total degree orthogonal expansion having order depending on s, we now provide a comparison 
in both of these cases. 

Comparison in the affine and non-affne polynomial cases, e.g., Examples 2.1 and 2.2. The terms (SG.2) 
from (47) and (SC.2) from (62) are asymptotic estimates of the number of coupled and decoupled finite 
element systems that must be solved by the SGFEM and SCFEM to construct the stochastic approxima¬ 
tion, respectively. For the coefficients from Examples 2.1 and 2.2, (SG.2) from (47) for the SGFEM is 
C9([log(l/er)] 7V ). Hence, our analysis shows that in these cases, the number of coupled finite element sys¬ 
tems in the SGFEM matrix is smaller than the number of decoupled finite element systems required by the 
SCFEM by a factor of (loglog(l/e)) Ar ~ 1 . This difference is enough to suggest that, if the condition numbers 
of both the preconditioned coupled system from the SGFEM and the preconditioned individual systems from 
the SCFEM are of the same order, then the SGFEM will outperform the SCFEM in terms of minimum work 
required to obtain a fully discrete approximation. In fact, whenever a[x,y) is a general non-affine, poly¬ 
nomial coefficient, e.g., Example 2.2, having fixed order r < 00, the complexity of matrix-vector products 
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involving K from (16) is approximately 0{ JhM p M r min{2 r , M^ r /2] }) from Remark 4.4. Hence, our analysis 
shows that when (D(M r min{ 2 r , }) < C , ((loglog(l/e)) Ar_1 ), which, in limit as e —> 0 , is always the 

case, and the condition numbers of both systems are of the same order, the SGFEM will outperform the 
SCFEM. However, in practical applications, it may require unrealistically small tolerance e to see this when 
r is large. 

Comparison in the non-affine, transcendental case, e.g., Example 2.3. In the case that a(x,y) is a 
non-afHne, transcendental functon of the random parameters, e.g., Example 2.3, the estimate for (SG.2) is 
0([log(l/e)] 3Ar ). Here our analysis shows that the number of coupled finite element matrices present in the 
SG system K r from (22) dominates the number of decoupled finite element systems needed by the SCFEM 
by a factor of [log(l/e)] 2iv . Note that, as in the case of the SGFEM, the term (SC.3) has a dependence on 
the condition numbers of the preconditioned finite element systems through the bound R = s\xp yer n(y). 
For the unpreconditioned systems A(y), the condition numbers can be bounded by /c(A (y)) < ( C K /h ) 2 
for every y € T, following from assumption (Al) and the quasi-uniformity of the mesh 7),. However, if we 
use the exact inverse of Po when preconditioning the SCFEM systems, the condition numbers are bounded 
independent of h and y £ T, since in this case R is independent of mesh size h and level L. Hence the work 
required by the PCG method when solving the SCFEM systems is dependent on £ only through the term 
log(l /s'). On the other hand, if we use the exact inverse of P when preconditioning the SGFEM system, the 
condition number h r can be bounded by 


h r < 


1 + T r 
1 — T ’ 


where r r and T r are defined in (49), hence depend on £ when r is chosen to satisfy r > r m j n from Lemma 
4.8. Figure 5 plots the condition numbers of both the unpreconditioned matrix K r and the preconditioned 
matrix P _1 K r with decreasing finite element mesh parameter h for the coefficient a(x, y) given in (67) from 
§6.3 with N = 4, L c = 1/2, and letting r = p with p increasing. There we see that the dependence on h 
has been removed by applying P _1 , but as p increases, we see a corresponding increase in the condition 
number. Other preconditioners than P may be used to reduce the dependence on r, e.g. [34], but then their 
associated costs must be accounted for in the work estimate (28) as well. However, even if the condition 
numbers of both the preconditioned coupled SG system and the preconditioned decoupled SC systems are 
of the same order, the additional work required to solve the coupled systems induced by the nonlinearity of 
the coefficient makes it difficult to see how the SGFEM can compete with the SCFEM. 

6. Numerical examples 

In this section, we provide illustrative numerical examples comparing the complexity of the SGFEM in 
the three cases of Examples 2.1, 2.2, and 2.3. We then compare these results with SCFEM and the results 
of the theoretical complexity comparison of the previous section. We solve the model problem (1), on the 
unit square D = [0, l] 2 . For a general coefficient a(x,y) we do not know the exact solution to (1). Hence 
we check the convergence against a “highly enriched” approximation, which we consider close enough to 
the exact one. To construct this “exact” solution u ex (x,y), we make use of the isotropic SCFEM based on 
Clenshaw-Curtis abscissas using the level L ex . We approximate the computational error for the SGFEM 
with orders p = 0,1,2,..., p rnax and for the SCFEM with levels L = 0,1,2,..., L max as 


||IE[£sg]IU°° « ||E[u ex - Uh, P ]\\e°° and ||E[£ S c]|U«> « l!E[u ex - Uft,L]||*> 


(63) 


where Uh, P and Uh,L are the fully discrete approximations (13) and (53), respectively, found by the PCG 
method, described in §3 and §5. In §6.3, we measure ||E[£sg] IU°° ~ ||E[u ex — u r hp \ ||^oo where u r h p denotes 
the solution of ( 22 ) with the projected coefficient a r (x,y). 

As stated in §3.4 and §5.2, we use PCG with the mean-based preconditioner for SGFEM and the level- 
zero preconditioner for the SCFEM. Hence, we believe this puts both methods at a similar starting point for 
comparison, if not providing a slight advantage for the SGFEM. With these choices, the complexity results 
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Figure 5: Condition numbers of both the unpreconditioned matrix K r and the preconditioned matrix P 1 K, with decreasing 
finite element mesh parameter h and r = p for the coefficient (67) from §6.3 with N = 4 and L c = 1/2. 


are presented in terms of the work estimates (28) and (58), respectively. The amount of work to reach a 
given error in PCG is also dependent on the tolerance used by the solver. If the tolerance is too small, we 
may see that the PCG method “over-resolves” the solution. To ensure that we do not over-resolve either 
solution, we set the tolerance of the solvers to be ||E[£sg]||^“/10 and ||E[£sc]||*°°/10 respectively, where 
these quantities are first estimated for each order p and level L using a tolerance of 1.0 x 10 -12 . In practice, 
we find that this does not affect the convergence results much. 

In all three examples, we use the SG approximation constructed in terms of the orthonormal Legendre 
polynomials {TpjpgA for given index sets A p . In the presentation of the results that follow, we use the 
following abbreviations. For the SGFEM, we use: “SG-TD” to denote the approximation in the total degree 
subspace 'P A TD (r) with A™ given in (12), and “SG-SM” to denote the approximation in the sparse Smolyak 
subspace Vasm (T) with A® M given in (12). For the SCFEM, we use: “SC-GL” and “SC-LJ” to denote the 
Smolyak approximation constructed on Gauss-Legendre abscissas and the Leja approximation constructed 
on Clenshaw-Curtis abscissas, both defined in terms of $td and m given in (51) and (55), respectively, and 
“SC-CC” to denote the Smolyak approximation constructed on Clenshaw-Curtis abscissas with g SM and m 
given in (51) and (54). 

6.1. Piecewise affine coefficients 

One common example in engineering and the physical sciences is that of isotropic thermal diffusion 
problem with a stochastic conductivity coefficient. Consider a partitioning of D = [0, l] 2 into 8 circular 
inclusions arrayed about 1 square inclusion as in Figure 6. We present the following example from [3], where 
the coefficient was given by 


8 

a(x,y) = b 0 (x) + ^2,VnXn{x), (64) 

n—1 

with b 0 = 1, and y n ~ U(— 0.99, —0.2). Here, \n are indicator functions corresponding to the 8 circular 
inclusions of radius r = 0.13. In this example, we also set the forcing term to be 

f[x) = 100xf(z), (65) 
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where F = [0.4, 0.6] 2 , is the square inclusion centered in D with side length 0.2. Figure 6 shows the expected 
value of the solution to this problem. To solve (1) with the coefficient (64) and forcing function (65), we use 
a piecewise linear finite element basis in the deterministic space over a nonuniform mesh Th- Here, the nodes 
of Th are adapted to the geometry of our problem, that is, we fix the nodes that lie on the boundaries of the 
inclusions in our domain. From this fixed boundary data, we then use the distmesh MATLAB program [28] 
to generate a non-degenerate triangulation that adequately resolves the details of our subdomain geometry. 
We further specify the subsets of the total set of nodes that belong to each geometric inclusion, and to the 
boundaries of the inclusion, so that the interface conditions for the coefficient may be correctly applied. The 
final mesh consists of 10,604 elements, 5,377 total nodes, and 5,229 unknowns. 



Figure 6: Left: a triangulation of the domain D with circular and square inclusions. Red nodes highlight the boundary of 
an inclusion or the domain D, blue nodes highlight nodes on the interior of an inclusion. Right: the expected value of the 
solution of (1) with stochastic conductivity coefficient (64). 


The coefficient (64) is an example of a coefficient a(x, y) having affine dependence on the parameters, e.g., 
Example 2.1. Figure 7 displays the convergence of the stochastic Galerkin and collocation methods against 
the total number of SDOF. For the SGFEM we take the SDOF to be the cardinality of the set A p used 
in constructing the fully discrete approximation Uh lP from (13) by solving (14), and for the SC method we 
take the SDOF to be the number of points ftG™' 9 corresponding to an index set A™' 9 used in constructing 
the fully discrete approximation Uh,L from (53). From the discussion of §5.3, we expect to see that the 
approximation obtained with the SGFEM requires fewer SDOF than the SCFEM to achieve the same error, 
and this is indeed the observed result. For example, both the SG-TD and SC-LJ approximations require 
the same number of SDOF, but the error of the SC-LJ approximation is much higher. This, of course, is 
a consequence of the estimate (59), where the errors of the SC approximations are bounded above by their 
respective Lebesgue constants against the best-approximation error in the space 

Figure 7 also displays the convergence of both methods in terms of error versus the total computational 
cost of solving the system with the work estimates of (28) and (58), respectively. Here, we compute the 
error in ||E[£sg]||{°o and ||IE[esc]|U°° as given in (60) and measure the cost as the number of O(Jh) matrix 
vector products required by both methods which are explicitly counted as 

Ktt? * ( ^ + E nnz(G r ) ) = N * (M p + M(p, r)) 

V reA r / 

in the code for the SGFEM and 2 * Y^l=i ^ter' i n the code for the SCFEM. Our analysis shows the work 
corresponding to the SG discretization for SG-TD is asymptotically bounded by O([log(l/e)] w ) while the 
analysis from [16] shows that the work corresponding to the SC discretization for SC-CC is asymptotically 
bounded by 0([log(l/e)] Ar [loglog(l/£)] Ar_1 ). This closely matches the results of the numerical experiments 
in Figure 7, where it can be seen that for polynomial order p > 2, the SG-TD approximation yields the best 
results with the least computational cost for this problem. 
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Figure 7: Left: Error versus SDOF in solving problem (1) with coefficient (64) and forcing (65). Right: Error versus 
computational cost with the work estimates given in (28) and (58) based on total number of matrix-vector products used by 
the CG method. 


6.2. Polynomial coefficients 

The next example we present is that of a polynomial function of the random paramters y 1 e.g. the 
coefficient from Example 2.2. We consider the following function 


a{x 1 y) = 5+ Y e lblrl ^r(x)y r , <; r (x) 

|r|<r 


sin (\r\-nxi) cos (\r\-iTX 2 ) if |r| is even, 
cos (|r|7ra;i) sin [\r\itX2) if |r| is odd, 


( 66 ) 


with y n ~ IA{— 1,1) for all n = 1,. .., TV and forcing term /( x) = 1 \/x £ D. For the results that follow we 
fix N = 4 and study the convergence of the SGFEM and SCFEM in the cases r = 1, 3, 7 in (66). As in §6.1, 
we set the finite element space for the spatial discretization to be the span of piecewise linear polynomials, 
but here we use a uniform triangulation of D with 4,934 elements and 2,340 spatial unknowns. 

Figure 8 displays the convergence of the SGFEM and SCFEM in terms of error versus the total com¬ 
putational cost of solving the system with the work estimates of (28) and (58). Here, we compute the 
error in ||E[esg]||.« 00 and IIEfesc]^ 00 as given in (63). As we increase the order r in (66), we see that 
the work for the SGFEM increases, corresponding to the decreasing sparsity of the matrix K from (16). 
Here, the work of matrix-vector multiplications with K are of the order 0(JhM p M ¥ min{2 r , M\ r /' 2 \}), where 
M P min{2 r , M^/i]} is a large constant that grows rapidly with r. As a result, we see that for r = 1, the 
SGFEM outperforms the other methods for p > 4. However, for f = 3, 7, the extra work of the matrix- 
vector multiplications of the coupled SG system dominates the overall convergence. We also observe that 
the convergence rate of the SGFEM does not change in these cases, as discussed in the comparison in §5.3. 

6.3. Transcendental coefficients 

The next example we present is that of a random coefficient defined in terms of the truncated Karhunen- 
Loeve expansion of the function log(a(cc, y) — a m i„), for a m i n > 0. This example represents a commonly used 
transcendental function of the physical and random parameters, e.g., Example 2.3, and is often presented 
in the context of enforcing the positivity of a(x,y) required in assumption (Al). Coefficients of this type 
are commonly found in groundwater flow models. For these models, the permeability can exhibit large 
variance within each layer of sediment, and as a result are better represented on a logarithmic scale. We 
recall the problem of solving (1) with a coefficient having one-dimensional (layered) spatial dependence and 
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cost 
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Figure 8: Error versus cost for solving problem (1) with coefficient (66) having r = 1 (left), r = 3 (middle), and r = 7 (right), 
with forcing f(x) = 1. The cost, given in (28) and (58), is based on total number of matrix-vector products used by the PCG 
method. 


a deterministic load /( Xi,X 2 ,u>) = 2 cos (a; 1 ) sin (# 2 ) as studied in [25, 26], where a(x,y ) was given by 


log(a(x, w) - 0.5) = 1 + yi(w) 


/ /— T \ N 

J + ^2,Cn<Pn{x)y n (uj) 


(67) 


<» - - (-<• f «">*<*>-{ 


Here, {yn{^)}^ = \ are independent random variables uniformly distributed in [—-\/3, -\/3] with zero mean and 
unit variance. For x\ £ [0,6], let L c be a desired physical correlation length for the random field a{x,y), 
chosen so that the random variables a(x 1 , ui) and a(x [, ui) become essentially uncorrelated for \x± — x\ | L c . 
Also, let L p = max{6, 2L C } and L = L c /L p . Expression (67) represents a possible truncation of a one- 
dimensional random field with stationary covariance, 

cov[log(a- 0.5)](xi, 2 : 2 ) = exp ^ j . 

Direct integration with the coefficient a(x,y ) from (67) yields a fully-block dense linear system K from (16) 
that is computationally infeasible to solve [14, 23, 34, 35]. The purpose of this example is to highlight the 
difficulties of obtaining a fully discrete approximation with the SGFEM in this case. 

As in the previous example in §6.2, we set the finite element space for the spatial discretization to be the 
span of piecewise linear polynomials and use a uniform triangulation of D with 4, 934 elements and 2, 340 
spatial unknowns. For the results that follow, we fix the truncation length N = 9 and correlation length 
L c = 1/64 in (67). To maintain sparsity of the SG system, we use the strategy of projecting the coefficient 
a(x,y) from (67) onto the space 7 5 A r (r), as in (18), where A r = A™ for the SG-TD approximation, and 
A r = A® M for the SG-SM approximation, obtaining the matrix K r from (22). We then increase r while p is 
fixed until the error in the solution stagnates, in practice finding that, for this problem, r = p is sufficient 
to guarantee the error of the projection does not exceed that of the solution, while maintaining sparsity of 
the linear system. 

Figure 9 compares the error versus SDOFs. There we see that for order p > 3, the SG-TD approximation 
provides the best approximation with respect to SDOFs. As discussed in §5.3, this is to be expected since the 
computational complexity of solving the coupled and decoupled systems is not taken into account. Figure 
9 also displays the convergence in error versus the total computational cost of solving the system with the 
work estimates of (28) and (58). Here however, the results show that the SGFEM requires significantly 
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SC-CC 

Level 

SC-CC 

Error 

Mat-vec cost 
of SC-CC 

SG-TD 

Order 

SG-TD 

Error 

Mat-vec cost 
of SG-TD 

0 

1.3626 x 10" 4 

2 

0 

1.3626 x 10” 4 

4 

1 

2.8884 x 10" 6 

218 

1 

3.9444 x 10~ 5 

152 

2 

6.3652 x 10" 8 

3,398 

2 

6.1427 x KT 7 

10,710 

3 

3.6021 x 10" 9 

28,638 

3 

2.8851 x 10~ 8 

213,010 

4 

1.4794 x 10" 10 

178,894 

4 

4.9210 x KT 10 

4,579,575 

5 

2.2869 x 10" 12 

944,220 

5 

8.9123 x KT 12 

49,089,051 


Table 1: Comparison of cost in matrix-vector products for solving problem (1) with coefficient (67) and forcing f{x i,X2,uj) = 
cos(xj) sin(iE 2 ) using the SC-CC and SG-TD approximations, with the strategy of picking the CG tolerance to be ||E[n ea; — 
u f h jjllfoo/10 for the SGFEM and ||E[u ra — /10 for the SC method. Cost in matrix-vector products for SG-TD method 

is given by (28) and for SC-CC is given by (58) normalized by the cost of a finite element matrix vector product. 


more work to obtain the same error than the SCFEM. We also observe the change in rate discussed in §5.3 
in this case, as the work required to solve (1) with the coefficient a(x, y) from (67) now depends on the 
order r of the projection used in the approximation of a(x, y). 

For the TD-SG approximation, this can be seen as a consequence of the fact that when r = p, the 
cost of solving (22) with the PCG method is of the order ©(J^M 8 , 2 -| -^2)> growing much more rapidly 
than the cost in the affine and polynomial coefficient cases, e.g., Examples 2.1 and 2.2, as we increase the 
order p. Table 1 shows the amount of work required to achieve an error on the order of 10“for some 
values of k in terms of the total number of matrix-vector products required by both the SC-CC and SG-TD 
approximations. 



Figure 9: Left: Error versus SDOFs in solving problem (1) with coefficient (67) and forcing f(xi,X2,w) = cos(.Xi) sin(.T 2 ). 
Right: Error versus cost with the work estimates given in (28) and (58). 


7. Conclusions 

In this work, we presented explicit cost bounds for applying the SGFEM to the solution of an elliptic 
PDE having both affine and non-affine random coefficients. To this end, we have conducted a rigorous 
counting argument for the sparsity of the linear system that results from the SG discretization with a global 
orthogonal basis defined on an isotropic total degree index set. Our analysis shows that when the coefficient 
is an affine or non-affine function of the random variables having fixed polynomial order, the computational 
cost of solving the coupled SG system grows linearly with the dimension of the polynomial subspace. In 
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these cases, the results only differ by a constant depending on the polynomial order of the random coefficient 
and the dimension of the parameter domain. 

On the other hand, when the coefficient is a non-affine, transcendental function of the random variables 
requiring an additional orthogonal expansion, our analysis shows that the computational complexity, no 
longer grows linearly with the polynomial subspace dimension. For such coefficients, we are able to provide 
bounds on the complexity that depend on the truncation order of the coefficient. These estimates imply 
that a truncation of the expansion should be used, when possible, though attention must be paid to the 
well-posedness of the resulting PDE. 

The analysis conducted herein motivates the study of the total computational complexity of obtaining 
fully discrete approximations with such methods. We have seen that, despite the fact that the SG method 
yields an approximation that is optimal in the L 2 sense for a given polynomial subspace, the associated 
computational costs of obtaining SG approximations are not optimal for all problems. Moreover, we have 
observed, both through theoretical comparisons and numerical examples, that changing the underlying 
polynomial subspace and method used for obtaining the fully discrete approximation can often yield a 
solution that requires far less work to obtain, but has the same error. 
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8. Appendix 

Proof (of Corollary 4.3). When N = 1 we denote r = r E No and note that from Theorem 4.1, 

c(r e) = / r even > e = r / 2 
|y#S(r, £) otherwise, 

with S(r, t) = {s € No : s = £, s < r} = {£} for i<r and 0 otherwise. We distinguish in cases: 

• Case r = 2k, k G No, 

1. when 0 < r < p, 



= (1 + p — k) — k(3k — 2p — 1) 
= 1 + p — 4 k 2 + 2 kp + k 2 
= (p — 2k + 1)(2 k + 1) + k 2 
= (p - r + l)(r + 1) + k 2 . 
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2. when p + 1 < r < 2 p, we have A— < k < p, so 


5 ( G r) = Y 

i=\r/2\ 


1 + p — A _ (1 + p - k 

p — £ J \ p — k 


E ^ 


1+p- £ 


=(!+?-*)+ e i 1+ „-~ e e ) 

e=k +1 ' " ' 

= (1 + p - k) + (p - k)(p - k + 1) 
= (1 +p — k ) 2 . 


3. when r > 2 p, then k > p, so 

r 

nnz(G r ) = c(r, 

t =|>/21 


1+p - A _ fl+p-k 
p — i ) \ p — k 


E ^ 


1+p- £ 


since p — k < 0 and l>k=>p—£<p — k< 0. 

• Case r = 2k + 1, k £ No, 

1. when 0 < r < p, then \r/ 2] = [(2 k + l)/2] = |"fc + 1/2] = k + 1, so 

-«v>= t ( 1 + p r e e 

1= fr/2] V P 1 ' e=k +1 \ P 1 ■ 


= —(1 + k)(3k — 2 p) 

= —4/c + 2p - 4fc 2 + 2fcp + k 2 + k 
= -2k(2k + 2) + p(2k + 2) + k 2 + k 
= (p — 2k) [2k + 2) + k 2 + k 
= (p — r + l)(r + 1) + k 2 + k. 


2. when p + 1 < r < 2p, then p/2 < k < p — 1/2, so 


<Gr>- £ «m>( 1 :!7')-*E ( , :'7‘)-»E 


t=+/21 


3. when r > 2p, then k > p — 1/2, so 


1 +p — £ 


1 +p — 


/G r ) = ^ c(r,J 

t=+/2l 


= (p - fc)(p — fc + 1). 


1+p : f ) = E=M)( 1+I, ;0-o, 

/ < 5 + v f-i l 


since k > p — 1/2 =+ p — £ < p — (k + 1) = p — k — 1 < p — (p— 1/2) — 1 = —1/2. 
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